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Abstract 

The purpose of this paper is to construct universal, auto-adaptive, localized, linear, polynomial 
(-valued) operators based on scattered data on the (hyper-)sphere S q (q > 2). The approxima- 
tion and localization properties of our operators are studied theoretically in deterministic as well as 
probabilistic settings. Numerical experiments are presented to demonstrate their superiority over 
traditional least squares and discrete Fourier projection polynomial approximations. An essential 
ingredient in our construction is the construction of quadrature formulas based on scattered data, 
exact for integrating spherical polynomials of (moderately) high degree. Our formulas are based on 
scattered sites; i.e., in contrast to such well known formulas as Driscoll-Healy formulas, we need not 
choose the location of the sites in any particular manner. While the previous attempts to construct 
such formulas have yielded formulas exact for spherical polynomials of degree at most 18, we are able 
to construct formulas exact for spherical polynomials of degree 178. 

Keywords: Quadrature formulas, localized kernels, polynomial quasi-interpolation, learning theory on 
the sphere. 

AMS classification: 65D32, 41A10, 41A25 



1 Introduction 



The problem of approximation of functions on the sphere arises in almost all applications involving 
modeling of data collected on the surface of the earth. More recent applications such as manifold matching 
and neural networks lead to the approximation of functions on the unit sphere S q embedded in the 
Euclidean space R 9+1 for integers q > 3 as well. Various applications in learning theory, meteorology, 
cosmology, and geophysics require analysis of scattered data collected on the sphere [3 IHl E] • This means 
that the data is of the form {(£, /(£))} for some unknown function / : S q — > R, where one has no control 
on the choice of the sites £. 

There are many methods to model such data: spherical splines, radial basis functions (called zonal 
function networks in this context), etc. However, the most traditional method is to approximate by 

"The research of this author was supported by Australian Research Council under its Centres of Excellence Program. 
tThe research of this author was supported, in part, by grant DMS-0605209 from the National Science Foundation and 
grant W911NF-04- 1-0339 from the U.S. Army Research Office. 
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spherical polynomials; i.e., restrictions of algebraic polynomials in q + 1 variables to S 9 . Apart from 
tradition, some important advantages of polynomials are that they are eigenfunctions of many pseudo- 
differential operators which arise in practical applications, and that they are infinitely smooth. Unlike 
in the case of spline approximation with a given degree of the piecewise component polynomials, global 
polynomial approximation does not exhibit a saturation property [5] Section 2, Chapter 11]; i.e., for 
an arbitrary sequence 6 n j 0, it is possible to find a continuous function on the sphere, not itself a 
polynomial, which can be approximated by spherical polynomials of degree at most n uniformly within 
S n , n > 1. In |20l 119) . we have shown how a good polynomial approximation yields also a good zonal 
function network approximation. In [22], we have shown that the approximation spaces determined by 
zonal function network approximation are the same as those determined by polynomial approximations. 

To illustrate the issues to be discussed in this paper, we consider an example in the case q = 1, 
or equivalcntly, the case of 27r-periodic functions on the real line. In this discussion only, let f(x) = 
Icosa;! 1 / 4 , x S R. In Figure Deleft), we show the log-plot of the absolute errors between / and its 
(trigonometric) Fourier projection of order 31, where the Fourier coefficients are estimated by a 128 point 
DFT. In Figure QJright), we show a similar log-plot where the Fourier projection is replaced by a suitable 
summability operator (described more precisely in (|3.1[) ). yielding again a trigonometric polynomial of 
order 31. It is clear that our summability operator is far more localized than the Fourier projection; i.e., 
the error in approximation decreases more rapidly as one goes away from the singularities at 7r/2 and 
37r/2. The maximum error on [37r/4, 57r/4] is 0.0103 for the projection, 0.0028 for our operator. Out 
of the 2048 points considered for the test, the error by the summability operator is less than 10~ 3 at 
38.96% points, the corresponding percentage for the projection is only 4.88%. In contrast to free-knot 
spline approximation, our summability operator is universal; i.e., its construction (convolution with a 
kernel) does not require any a priori knowledge about the location of singularities of the target function. 
It yields a single, globally defined trigonometric polynomial, computed using global data. Nevertheless, 
it is auto-adaptive, in the sense that the error in approximation on different subintervals adjusts itself 
according to the smoothness of the target function on these subintervals. In [23] [24] . we have given a 
very detailed analysis of the approximation properties of these operators in the case q = 1. 
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Figure 1: The log-plot of the absolute error between the function x h-> | cosxl 1 / 4 , and (left) its Fourier 
projection (right) trigonometric polynomial obtained by our summability operator, where the Fourier 
coefficients are estimated by 128 point FFT. The order of the trigonometric polynomials is 31 in each 
case. The numbers on the x axis are in multiples of 7r, the actual absolute errors are 10 v . 

Our computation based on a 128 point DFT implies that the values of the function are available at 
128 equidistant points. If only a scattered data is available, the following method is often used (especially 
in the context of approximation on the sphere) to estimate the values needed for the DFT. For each 
point £, we consider the nearest point of the form 2-7rfc/128, and imagine that the value of / at this point 
is /(£), taking averages in the case of multiplicities, and interpolating in the case of gaps. If we use 
our summability operator, estimating the Fourier coefficients in this way, then the maximum error on 
[37r/4, 57r/4] is 0.0357, and the proportion of points where the error is less than 10 -3 is 7.08%. It is clear 
that a careful construction of quadrature formulas is essential to obtain good approximation results. 

The purpose of this paper is to construct universal, auto-adaptive, localized, linear, polynomial (- 



2 



valued) operators based on scattered data on § 9 (q > 2) and to analyse their approximation properties. An 
essential ingredient in our construction is the construction of quadrature formulas based on scattered data, 
exact for integrating spherical polynomials of (moderately) high degree, and satisfying certain technical 
conditions known as the Marcinkiewicz-Zygmund (M-Z) conditions. Our construction is different from 
the usual construction of quadrature formulas (designs) studied in numerical analysis, where one has a 
choice of the placement of nodes. In [5T|, we had proved the existence of such quadrature formulas for 
scattered data. However, previous efforts to compute such formulas did not yield exactness beyond degree 
18 polynomials. This was a severe limitation on the practical applications of our theoretical constructions. 
We will show that a very simple idea of solving a system of equations involving a Gram matrix yields 
surprisingly good results, in particular, quadrature formulas exact for integrating polynomials of degree 
as high as 178. Gram matrices are typically ill-conditioned. However, we will show both theoretically 
and numerically that the ones which we use are, in fact, very well conditioned. We will introduce another 
algorithm of theoretical interest to compute data dependent orthogonal polynomials, and use these to 
compute the quadrature formulas in a memory efficient manner. To the best of our knowledge, this is 
the first effort to extend the univariate constructions in Gautschi's book [11] to a multivariate setting. 
Considering that computation of classical spherical harmonics is a very delicate task, requiring many 
tricks based on the special function properties of these polynomials for a stable computation, it is not 
expected that our computation of data dependent orthogonal polynomials with no such special function 
properties would be stable. In describing this algorithm, we hope to stimulate further research in this 
interesting direction. We note that even if this algorithm is not as stable for high degrees as the other 
algorithm, it yields satisfactory quadrature formulas exact for integrating polynomials of degree 32. Most 
importantly, our new found ability to compute quadrature formulas for moderately high degrees allows 
us to offer our operators as a viable, practical method of approximation, even superior to the commonly 
used methods of least squares and Fourier projection as far as localized approximation is concerned. 

An additional problem is when the available values of the target function are noisy. One may assume 
that the noise is an additive random variable with mean zero. It is also routine in learning theory to 
assume that the random variables have a bounded range. This assumption is usually satisfied with a high 
probability even if the random variables do not actually have a bounded range. However, one does not 
typically know the actual distribution of these random variables. We obtain probabilistic estimates in 
this setting on the global and local approximations by our operators. To underline the practical utility 
of our operators, we use them for modelling the MAGSAT data supplied to us by Dr. Thorsten Maier, 
obtaining results comparable to those obtained by other techniques. 

In Section[2l we review certain facts about spherical polynomials, the existence of quadrature formulas 
to integrate these, a few properties of the quadrature weights, and certain polynomial kernels which we will 
need throughout the paper. In Section[3j we study the approximation properties of the linear polynomial 
operators. The new results here are Theorems 13.11 and 13.21 The first parts of these theorems were proved 
essentially in |18j . but not stated in the form given here. In order to apply these operators in practice, one 
needs quadrature formulas exact for high degree spherical polynomials. Explicit algorithms to construct 
such formulas are described in Section |4j The new results in this section are Theorems 14.11 and 14.21 
Numerical results are presented in Section [5j and the proofs of all new results are given in Section [6l The 
paper is a result of a long process, involving discussions with a number of mathematicians. In particular, 
it is our pleasure to acknowledge the support and encouragement of Mahadevan Ganesh, Thorsten Maier, 
Volker Michel, Dominik Michel, Ian Sloan, and Joe Ward. We are also grateful to the two referees and 
Fred Hickerncll for their many useful suggestions for the improvement of the first draft of this paper. 



2 Background 



In this section, we review some known results regarding spherical polynomials and localized polynomial 
kernels. 



2.1 Spherical polynomials 

Let q > 1 be an integer, S 9 be the unit sphere embedded in the Euclidean space i.e., 

S« := {(si,... , Xq+i) € R 9+1 : — 1}' ana - be its Lebesgue surface measure, normalized so 

2jT (q+l)/2 

that fiq(§ q ) = 1. The surface area of § 9 is —77 ttt- F° r <5 > 0, a spherical cap with radius S and 

r((g+l)/2) 

center Xo € § 9 is defined by 

S^(x ) := {xeS' : arccos(x • x ) < 5}. 
If 1 < p < 00, and / : S 9 — > M is measurable, we write 



IP 



_ /{/ s J/(x)|^ 9 (x)} 1/p , if 1 < P < 
ess sup x6S , |/(x)|, if p = 00. 



00, 



The space of all Lebesgue measurable functions on S 9 such that ||/|| p < 00 will be denoted by L p , with 
the usual convention that two functions are considered equal as elements of this space if they are equal 
almost everywhere. The symbol C(S 9 ) denotes the class of all continuous, real valued functions on S 9 , 
equipped with the norm || o W^,. 

For a real number x > 0, let II 9 denote the class of all spherical polynomials of degree at most x. 
(This is the same as the class II 9 , where n is the largest integer not exceeding x. However, our extension 
of the notation allows us, for example, to use the simpler notation n 9 y 2 rather than the more cumbersome 
notation n 9 n ^ 2 j .) For a fixed integer £ > 0, the restriction to § 9 of a homogeneous harmonic polynomial 
of exact degree £ is called a spherical harmonic of degree £. Most of the following information is based 
on [55], [551 Section IV. 2], and [5J Chapter XI], although we use a different notation. The class of all 
spherical harmonics of degree £ will be denoted by H 9 . The spaces H 9 are mutually orthogonal relative 
to the inner product of L 2 . For any integer n > 0, we have IT 9 = (J)" =0 H 9 . The dimension of H 9 is 
given by 

, M + g -l/< + g -l\ 
d\ :=dimH« = { t+q-1 V I ) (2.1) 

1, if £ = 0. 

and that of II 9 is Y^e=o^£ = Furthermore, L 2 — i 2 -closure{ ©£^ -^?}- Hence, if we choose an 

orthonormal basis {Y^fc ■ k — 1, . . . , df } for each H 9 , then the set \Yi t k '■ £ = 0, 1, . . . and k = 1, . . . , d\} 
is a complete orthonormal basis for L 2 . One has the well-known addition formula [25] and [5J Chapter 
XI, Theorem 4]: 

d " 1 91- 1 r(n/9\ 2 

£y Afe (x)Y,, fe (C) = i^iL w(1)w(x . c)j £ = 0,1,..., (2.2) 

where pi := pf^" 1 1 ' q ^ 2 1 ' is the orthonormalized Jacobi polynomial with positive leading coefficient: 



Pe(t) Pk (t)(l-t 2 ) q / 2 - 1 dt = 



1, if£ = fc, 
0, otherwise. 



In particular, for x 6 § 9 , ^ = 0, 1, • • •, 

d q d q 

X>A« = 2g ' lr y) 2 p£(1) 2 = f J2Y 2 k (CW q (C) = d\. (2.3) 

fc=l 1 W ^ fc=l 

2.2 Localized polynomial kernels 

Let h : [0, 00) — > K be a compactly supported function, and i > 0. We define for ueE, 

? g_1 rfo/2 N i 2 00 

Mh;u) := (2.4) 

^ £=0 
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and define $t(/i; u) = if t < 0. 

In the sequel, we adopt the following convention regarding constants. The letters c, c\, ■ ■ ■ will denote 
generic, positive constants depending only on the dimension q and other such fixed quantities in the 
discussion as the function h, the different norms involved in the formula, etc.. Their value will be different 
at different occurrences, even within the same formula. The symbol A ~ B will mean cA < B < c\A. 

The following proposition summarizes some of the important properties of the kernels defined in (|2.4|) . 



Proposition 2.1 Let S > q be an integer, h : [0, oo) — > M be a S times iterated integral of a function of 
bounded variation, h(x) — 1 for x £ [0, 1/2], h(x) — for x > 1, and h be non-increasing. Let x£§'. 
We have for every integer n > 0, Q n (h; o • x) 6 Tl^, and 

y$ n (fc;x.0P(OdMO=-P(x), p en n/2 . (2.5) 



Further, 



sup / £.014^,(0 =sup / |$„(/i;x-£)|d/ig(£) 

ra>l, C£S ? J «>1 J 

■sup / |# n (fe;u)|(l-u 2 )«/ 2 - 1 du<oo, (2.6) 



C6S 

2W 2 

r(«/2) 



/ 



r(g/2) 

~ n 9 ~ max |$„ f/i: x 

and for every ^ £ S', x 



n q ~ max |$ n (/i;x • 01 = 1)|. (2.7) 



Except for (|2.7p , all parts of Proposition ^. ll have been proved and verified repeatedly in [T7l[T8lfT2l[T9] . 
We will sketch a proof of this proposition, mainly to reconcile notations. 

PROOF of Proposition 12.11 The equation (|2.5p and the first two equations in (|2.6[) are clear. The last 
estimate in (12 .6[) follows from p~7[ Lemma 4.6] with following choice of the parameters there: a — (3 = 
q/2 — 1, h v = h(v /n), where we observe that by a repeated application of the mean value theorem, 

oo 

^(i/ + l) s |A r /i(j//ra)| < cn s - r+1 , set, r, n = 1, 2, • • • , 



u=0 



where A r is the r— th order forward difference applied with respect to v. Similarly, the estimate v _ 
follows from [I7[ Lemma 4.10] with same parameters as above, S in place of K in [T7], and y = x • £ (cf 
Appendix to [12]). We prove (|2.7p . The first equation is a consequence of the rotation invariance of \i q 
In view of the addition formula (12.21) . 



n 

$ n (h; x • = J2 Kt/n) J2 Yt, k (x)Y tik (£). (2.9) 

1=0 k=l 

It follows using ([23]) and the facts that h(t/n) = 1 for t < n/2, < h(t) < 1 for t e [0,oo), that 

„ n a i n 

/ $ n {h-, x • o 2 d^ Q (o = E mw 2 E y ^( x ) 2 = E m^) 2 4 - « 9 - 

£=0 fe=l fcO 
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Similarly, using Schwarz inequality, l|2.2[) . (|2.3|) . and the fact that h(£/n) > 0, 
$„(/i;l) = |$„(/j;x-x)| < sup |$„(/i;x-£)| 

n a £ a i n 

< J2 y e,k (x) 2 } 1/2 {E Y eA0 2 } 1/2 = E = !)■ 

£=0 fe=l fe=l ^=o 

Since <2| ~ < < 1, and h(£/n) = 1 for t < n/2, the above two estimates lead to (|2~7| . □ 

In the remainder of this paper, h will denote a fixed function satisfying the conditions of Proposi- 
tion o 



2.3 Quadrature formulas 

Let C be a finite set of distinct points on A quadrature formula based on C has the form Q(f) = 
S^gc w ^f{^)i where iff , £ S C, are real numbers. For integer n > 0, the formula is exact for degree n if 
Q(P) — L q Pdfj, q for all P € II*. It is not difficult to verify that if Q n {f) = J2 w £ n f(£n) is a sequence 
of quadrature formulas, with Q n being exact with degree n, then Q n (f) — > J fdfi q for every continuous 
function / on S 9 if and only if ^ \ w £ n \ —■ c , with c being independent of n. In the sequel, we will assume 
tacitly that C is one of the members of a nested sequence of finite subsets of S 9 , whose union is dense 
in S q . All the constants may depend upon the whole sequence, but not on any individual member of 
this sequence. Thus, a formula Q will be called a bounded variation formula if X^ee \ w t\ — c ' with the 
understanding that this is an abbreviation for the concept described above with a sequence of quadrature 
formulas. 

Definition 2.1 Let m > be an integer. The set C admits an M-Z quadrature of order m if there exist 
weights such that 

f P(x)dMx)=5>^(0, Fen '- ( 2 - 10 ) 

and 

\ i/p 

EKII^W <c||P|| p , PenL, i<p<». (2.11) 

The weights will be called M-Z weights of order m. The condition H2.11\) will be referred to as the 
M-Z condition. 

If C admits an M-Z quadrature of order m, and {w^} are the weights involved, it is clear from using 
(|2.11|) with the polynomial identically equal to 1 in place of P that J^^ec \ w c\ — c - Further, if £ G C, 
then applying (|2.1ip with p = 2 and Q m (h; ( ■ °) in place of P, we obtain for M-Z weights of order m: 

\w C \$ m {h] l) 2 < E \ w t\*™(h; C • 2 < c f $ m (h; C • x) 2 d Mg (x). 

The estimate (|2.7[) now implies that for all M-Z weights {w^} of order m, 

|iu s | < cm~ q , £eC. (2.12) 

In |21j . we proved that every finite set C C §' admits an M-Z quadrature with an order depending 
upon how dense the set C is. This density is measured in terms of the mesh norm. The mesh norm of C 
with respect to a subset K C § 9 is defined to be 

8 C (K) := sup dist(x,C). (2.13) 
The following theorem summarizes the quadrature formula given in [21j . 



G 



Theorem 2.1 There exists a constant ct q with the following property. Let C be a finite set of distinct 
points on Efl , and m be an integer with m < a q (Sc(S q ))^ 1 . Then C admits an M-Z quadrature of order 
m, and the set {w^} of M-Z weights may be chosen to satisfy 

|U : w ( ^ 0}| ~ m« ~ dim(nlJ. (2.14) 



3 Polynomial operators 

For t > 0, we define the summability operator a\ by the formula 

<(ft.;/,x)=/ f(0^t(h]^-Od^ q (0 = J2 h ^J2f^ k ^ k ^' / G L 1 , x € § ? . (3.1) 

•' Sq 1=0 k=l 

(It is convenient, and customary in approximation theory, to use the notation a^(h;f, x) rather than 
of (ft; /)(x).) Although we defined the operator for L 1 to underline the fact that it is a universal operator, 
we will be interested only in its restriction to C(S 9 ). If / : E> q — > R is a continuous function, the degree 
of approximation of / from 11^ is defined by 

E x {f)= inf l|/-P||oo. 
Pen' 

It is well known pHQI] that for all integer n > 1, and / € C{S q ), 

E n (f) < 11/ - <{h\ f)\\oo< cE n/2 (f). (3.2) 
Following [18] . we now define a discretized version of these operators. 

If C C is a finite set, W = {u^j^gc an d Z = {z^}^ 6 c are sets of real numbers, we define the 
polynomial operator 

W£Z£G>t(h;x.- i e K, x£§'. (3.3) 

If / : §« R, and ^ = /(^), ^ e C, we will write cr t (C, W; h; f, x) in place of cj f (C, W; ft,; Z, x). In [18], 
we had denoted these operators by o~t(v; h, /), where i/ is the measure that associates the mass with 
£ € C. In this paper, we prefer to use the slightly expanded notation as in (|3.3|) . If n > 1 is an integer, 
C C § ? is a finite set that admits an M-Z quadrature of order n, W is the set of the corresponding M-Z 
weights. Then it is shown in [18l Proposition 4.1] that 

E n (f)< ||/ -<t„(C,W; ft; /)|U < cB n/a (/), /eC(§«). (3.4) 

In this paper, we will be especially interested in the approximation of functions in the class W r , r > 0, 
comprised of functions / € C(S 9 ) for which E n (f) = 0(n~ r ), n > 1. A complete characterization of the 
classes W r in terms of such constructive properties of its members as the number of partial derivatives 
and their moduli of smoothness is well known [551 US]- In view of (|3.2[) . / S W r if and only if 

||/||w r := ll/lloo + su P 2"la 2 *„(ft; /) - ^(h; f)^ < oo. 

n>l 

In practical applications, the data is contaminated with noise. Therefore, we wish to examine the 
behavior of our operators based on a data of the form {(£, /(£) + e^)}, where are independent random 
variables with unknown probability distributions, each with mean 0. If the range of these random variables 
is not bounded, one can still assume that the probability of the variables going out of a sufficiently 
large interval is small. Hence, it is customary in learning theory to assume that the variables have 
a bounded range, so that one may use certain technical inequalities of probability theory, known as 
Bennett's inequalities; see the proof of Lemma RT21 below. 
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In the statements of the theorems below, we use three parameters. The symbol M denotes the number 
of points in the data set; we assume that the set admits an M-Z quadrature of order to, and the degree n 
of the polynomial approximant a n (C, W; h; f) is determined in terms of to. For theoretical considerations 
where one is not concerned about the actual numerical constructions of the quadrature weights, one may 
imagine a data set C with M := \C\ ~ 8c(§ q )~ 9 , and assume that the weights W are as guaranteed by 
Theorem l2.ll If so, then we may take M ~ m q ~ 6c{S q )~ q in the discussion in this section. For example, 
the estimates (|3.7p and (|3.8[) below can then be expressed in terms of the number of samples respectively 
as follows: 

\\f-a n (C,W;h;f)\\ 00 <cM- r/q , with some n ~ M 1 / 9 , (3.5) 

and 

Prob (\\a n (C, W; h; Z) - /IU > Cl ^^+1) ) ^ C 2 M ~ C > with some n ~ (M / \ogM) 1 ^ 2r+q \ 

(3.6) 

Theorem 3.1 Suppose that m > \ is an integer, C = {^j}jLi admits an M-Z quadrature of order m, 
and let W be the corresponding quadrature weights. Let r > 0, / G W r , ||/||w r = 1- 

(a) For integer n < m, we have 

\\f-tr n (C,W;h; <cn- r . (3.7) 

(b) For j — 1,---,M, let tj be independent random variables with mean and range [—1,1], Z = 
{tj + /(£;)}■ If A > and n > 1 is £/ie greatest integer with (A + q)n 2r+q logn < c^m q , n < m, then 

Prob (||a n (C,W;^;Z)-/|| 00 >cin- r ) <can~ A . (3.8) 

Here, the constants C\,C2, C3 are independent of the distribution of the variables ej . 

We now turn our attention to local approximation by our operators. In the sequel, if K C § 9 , 
f:K^R, then H/H^ x : — sup xg ^- |/(x)|. If xo G a function / is defined to be r-smooth at xo if 
there is a spherical cap §^(xo) such that ftp £ W r for every infinitely diffcrcntiablc function <f> supported 
on §*(xo). We have proved in [TBI Theorem 3.3] that / is r-smooth at a point Xo if and only if there is 
a cap §^(xo) such that 

\\o-* 2n (h;f) - ^-i(/i;/)||oo,s|(*o) = 0(2~ nr )- 

Accordingly, if if is a spherical cap, we may define the class W r (K) to consist of / € C(S 9 ), for which 

ll/llw r (K) := ||/||oo + sup2" r ||a5„(/ l ;/)- ( 7*„- 1 (/i;/)|| 0O ,K < 00. 
n>l 

Theorem 3.2 Suppose that m > 1 is an integer, C = {^j}jLi admits an M-Z quadrature of order m, 
and let W be the corresponding quadrature weights. Let < r < S — q, K' C K be concentric spherical 
caps, f g C(Si), 11/11 Wr{K) = 1. 

(a) For integer n, 1 < n < m, 

\\f - a n (C,W;h; m^K, < cn- r . (3.9) 

(b) For j = 1, • • • , M, let ej be independent random variables with mean and range contained in [—1, 1] , 
and Z = {e,j + /(£;)}■ If A > and n > 1 is the greatest integer with (A + q)n 2r+q logn < C3TO 9 , n < to, 
then 

Prob (K(C,W;^Z)-/|| 00|Jf / >cm- r ) < c 2 n" A . (3.10) 
Here, the constants c±,C2, C3 are independent of the distribution of the variables 6j . 
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4 Construction of quadrature formulas 



In this section, we describe two algorithms to obtain bounded variation quadrature formulas associated 
with a given finite set of points C C S 5 . Both of these constructions can be described in a very general 
setting. Since this also simplifies the notations and ideas considerably by avoiding the use of real and 
imaginary parts of a doubly indexed polynomial Yt t k, we will describe the algorithms in this generality. 

Let fl be a nonempty set, /i be a probability measure on f2, C C fl, yi, 1/2, • • • be a complete orthonormal 
basis for L 2 (fl, fi), where y% = 1, and Vk denote the span of j/i, • • • , yk- Let v be another measure on f2, 
and (o, o) denote the inner product of £ 2 (f2, v). For integer N > 1, the Gram matrix Gn is an N x N 
matrix, defined by (Gjv)^.fe = (ye,Vk) — (Gjv)fc^, 1 < k,£ < N. We wish to find a weight function W on 
f2 such that J n Pd[i = J n PWdv for all P e Vjv for an integer N for which Gat is positive definite. 

For the applications to the case of quadrature formulas for the sphere, Q = S q , fi = n q , and y^s are 
the orthogonal spherical harmonics, arranged in a sequence, so that y\ = 1, and all polynomials of lower 
degree are listed before those of a higher degree. To include all polynomials in IT 9 , we need N = <i 9+1 . 
There are many possibilities to define the measure v. The simplest is the measure v MC that associates 
the mass 1/\C\ with each point of C. A more sophisticated way to define the measure v is the following. 
We obtain a partition of S 9 into a dyadic triangulation such that each triangle contains at least one point 
of C. We choose only one point in each triangle, and hence, assume that each triangle contains exactly 
one point of C. We define the measure v TR to be the measure that associates with each £ G C the area 
of the triangle containing £. 

One of the simplest ideas to compute the quadrature weights is the following. Let N be an integer for 
which Gn is positive definite. If P — a jVjt then J n Pd/i — 01. Also, the vector a = (ai, . . . , &./v) T 

satisfies the matrix equation 

G A ra=((P, 2/1 ),...,(P, 2/7V )) T , 

so that 

N 

/ Pdn = ^(GaO^P,^) = (P^GaO^). (4.1) 

^ a fc=i fe 

In the setting of the sphere, this gives the following quadrature formula: 

f Pd N = j2 p (o (k{£})E( g *)i>^)) =: E w f Qp &- (4-2) 

,7s ' ! cec I fc=i J {eC 

We formulate this as 

Algorithm LSQ 

Input: The matrix Y = (yk(£,))> k = 1, - ■ ■ ,N (optional), and the vector v = ({£}). 

1. Solve Fdiag(v)y T b = (1,0,..., 0) T . 

2. Return Q = ^({^}) J2k=i hVkit). 

We observe that Gn = Ydi&g(\r)Y T . It is clear that the matrix Gat is always positive semi-definite; 
the assumption that it is positive definite is equivalent to the assumption that no clement of V/v vanishes 
identically on C. If C and {^({C})} satisfy the M-Z inequalities, Theorem l4.1l below shows that Gat is well 
conditioned. Assuming that the matrix Y is input, the time to compute Gat is 0(7V 2 |C|) and the space 
requirement is 0(N 2 ). (In the case of the sphere S 9 , we need N = d 9 " 1 " 1 = (D(n q ) to compute formulas 
exact for degree n.) The vector b in Step Q] can be found using such iterative methods as the conjugate 
residual method. We refer to [6j for a more detailed analysis of this method. Using this approach, 
the matrix Y and Gat need not be stored or precomputcd, but the product of the matrix Gat with an 
arbitrary residual vector r needs to be computed. This observation results in a substantial saving in the 
time and memory complexity of the algorithm when the results are desired only within a given accuracy. 
For the unit sphere § 2 , when N = (n + l) 2 , the product Gatt = Fdiag(v)F T r can be computed within 
an accuracy e using a recent algorithm of Keiner [T?] using 0(rt 2 (log?i) 2 + log(l/e)|C|) operations, where 
e is the accuracy of the method. 
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One way to interprete this algorithm is the following. Let / : E> q — > R, and P be the solution to the 
least square problem 

P = argmm{(/-Q,/-Q) : QeV N }. 

If f is the vector ({/,%•)), then P — Y^ji^t^^jVj- The quadrature formula with weights w^ S ® thus 
offers Jg, Pdfiq as the approximation to J sq fd\i q . The weights w^ S ® also satisfy a least square property 
among all the possible quadrature formulas, as shown in Lemma 16. If a). We summarize some of the 
properties of the weights w^ S ® in the following theorem. 

Theorem 4.1 Let n > 1 be an integer, N = df^ 1 , C be a finite set of points on E> q and v be a measure 
supported on C. Let v% := £ € C, and 




ci\\P\\ P <\J2 v i\ p (0\ p \ <C2\\P\\ P , Pen', i< P <oo, (4.3) 

(a) For the Gram matrix Gn , the lowest eigenvalue is > c\, and the largest eigenvalue is < c\, where c\, C2 
are the constants in |^.ff| ) with p — 2. In particular, Gn is positive definite. Moreover, X^ec l^f I — c - 

(b) // 

J P 2 dv - J P 2 dii q <^J P 2 dfx q , P S II* , (4-4) 

then \w^ SQ \ < cv^, £ G C. In particular, the weights {w^ S< ^} satisfy the M-Z condition. 

(c) Let M > 1 be an integer, C be a independent random sample of M points chosen from the distribution 
fj,q, and A, r\ > 0. Let = l/M, £ G C. There exists a constant c = c{A) such that if n > 2 is an integer 
with M > cn q log n/rj 2 , then 



Prob 



P^di/ - / P 2 d/x 



<■</ 



> 7? / p^m 3 , Pen' < Cl n- A . (4.5) 



In particular, if M > cn 3q logn 7 then the condition \4-.J$ is satisfied with probability exceeding 1 — c\n A . 

One disadvantage of the algorithm LSQ is that one needs to know the value of N in advance. We 
now describe an idea which has the potential to avoid this problem. In the case when 51 is a subset 
of a Euclidean space, and the t/j's are polynomials, with y\ denoting the constant polynomial, one can 
construct a system {tk} of orthonormalized polynomials with respect to v using recurrence relations. 
Recurrence relations for orthogonal polynomials in several variables have been discussed in detail by 
Dunkl and Xu [H Chapter 3]. In contrast to the viewpoint in [3], we may depend upon a specific 
enumeration, but require the recurrence relation to have a specific form described in Theorem 14. 21 below. 
This form allows us to generalize the ideas in Gautschi's book [11] Chapter 2] in our context. 

To describe our ideas in general, let fl C R 9+1 , U\,U2, - ■ ■ be the lexicographic enumeration of the 
monomials in q + 1 variables, so that u\ is the monomial identically equal to 1, the restrictions of itfc's 
to Q are linearly independent, and Vk = span {u\, • • • , Uk}- It is not difficult to see that for every integer 
k > 1, there is a minimal index p(k) such that there exists a monomial fk of degree 1 with 

/fcWp(fc) = Uk+i, k = l,2,---. (4.6) 

We now let, for each k = 1, 2, • • •, {y±, • • • , yk} be a basis for Vk orthonormal with respect to /i, N > 1 be 
an integer for which the Gram matrix Gn is positive definite, and for each k = 1, • • • , JV, {t\, • • • , tk} be 
a basis for Vk orthonormal with respect to v. Clearly, any polynomial P £ Vn can be written in the form 



P(x)= f p(C)5>(x)t fe (c>MC), 



and consequently, one gets the "quadrature formula" 



J p(x)d M (x) = J p(o I (/ **w^(x)) t k (o I du(o. 



(4.7) 
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In this discussion only, let tk ='■ Y]j Ck,jUj, and the matrix (ckj) be denoted by C. The condition that 
ti, • ■ • ,t n is an orthonormal system with respect to v is equivalent to the condition that CGnC t = I, 
where I is the N x N identity matrix. Hence, G^ 1 = C T C. Moreover, j t k d^ = c k ,i for k = 1, • • • , N, 
and hence, we conclude that 

Y ( / ifc(x)d/i(x) J t k =^2^2c k ,iCkjyj = Y(G N )ijy 3 - 

k VJ / 3 k 3 



Thus, the quadrature weights in (|4.7[) are the same as those in (14. ip . 

First, we summarize the various recurrence relations in Theorem 14.21 below, although we will not use 
all of them. We will denote the (total) degree of u k by Dk, and observe that Dk is also the degree of yk 
and tk, Dj < j, and D p( j^ = D k +\ - 1. 

Theorem 4.2 There exist real numbers Skj, fk,j, Ak > 0, and a linear polynomial f k , such that 

fkVp(k) = Vk+i - ^ rkjVj, fktp(k) = Aktk+i - ^ s k,j t 3- ( 4 - 8 ) 

O k + 1 -2<D j <D k D k + 1 -2<D 3 <D k 

More generally, if P is any linear polynomial, there exist real numbers rkj(P) such that 

PVk= Y, r *AP)V3i (4-9) 

D k -l<D 3 <D k + l 

We have tk = J2j c k,jVj> where 

AkCk+1,1 = < Y r Lm(fk)Cp(k), m + Y S k,3 c 3,£ ( ■ ( 4 -10) 

(De.-l<D m <D e +l D k+1 -2<D 3 <D k J 

In the context of the sphere § 9 , we will compute f^'s using (|4. 8[) . and compute J tkdfi q using a known 
quadrature formula. The resulting algorithm, Algorithm REC, in the context of the sphere is summarized 
below. This algorithm is similar to the Stieltjes method in Gautschi's book [TlJ Section 2.2]. Even though 
it is feasible to carry out the algorithm for as large an N as the data allows, and to find this value of 
N during run time, it is still desirable from the point of view of numerical stability to limit the largest 
iV from the outset. Accordingly, in describing the following algorithm, we stipulate that the quadrature 
formula is to be computed to be exact only for polynomials in Vn for the largest possible N < L for 
some integer L > 1. We assume further that we know another quadrature formula (for example, the 
DriscolHIealy formula [3]) exact for polynomials in Vl- 

kP(0 = f Pd^ qi P G Vl- (4.11) 

Algorithm REC 

Input: An integer L, the sequence p(k), k = 1, • • • ,L, sets C, C*, weights (A^)^ e c» so that (|4.11[) holds, 
the values {yj(£)}seC, {%(C)}cec* for j = 1,2,3,4, and the values f k (£), f k ((), k = l,---,L. 

1. Using Gram-Schmidt procedure, initialize t\, ti, ts, ti, both for points in C and in C* , and initialize 
N = 4. 

2. For k = 1, • • • , 4, let lk = E C ec* ^MO- 

3. For each £ 6 C, initialize — X)t=i 7fctfc(0- 

4. For k = 4, 5, • • • (so that the degrees are at least for all polynomials entering in the recursions) 
and while N < L, repeat steps [5H5] below. 
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5. For j with D k +i - 2 < Dj < D k , set 

Sk,j = (fkt p (k),tj). 

6. Define Tk+\ by 

Tk+i = fktp(k) - E s fe.i*J- 

r> fc + 1 -2<_D J <D fc 

both for points in C and points in C*. if Ik+i = (Tk+i,Tk+i) = 0, then stop, and set N = k. 

1/2 

Otherwise, define t^+i = T k +i/I k ' +1 . 

7. Set 7 fc+ i = Ecec* A C f fc+i(0- 

8. For each ^eC,w i =w £ + - /k+1 t k+1 (^) 1 k = k + l, N = N+1. 

In the case of the sphere § ? , we take L = for some integer n > 1. The number of j's with 

ffe+i -2<Dj<D k ,l<j,k<L is ©(n^ 1 ). In this discussion only, let M = \C\ + \C*\. Consequently, 
Steps [5] and [6] require 0(Mn q ~ 1 ) operations. Since the remaining two steps in the loop take O(M) 
operations, the loop starting at Step 4 require 0(Mn 2q ~ 1 ) operations. Finally, we observe that in 
implementing the above algorithm, one need not keep the whole matrix only the rows corresponding 

to three degrees are required in any step. In particular, the memory requirement of this algorithm is 
OfMA'" 1 ). 



5 Numerical experiments 

The objective of this section is to demonstrate and supplement the theoretical results presented in Sec- 
tions [3] and |4j 

Our first set of experiments illustrates the algorithms LSQ and REC. The experiments were con- 
ducted over a long period of time, many of them long before we started to write the paper. Hence, the 
normalizations for the spherical polynomials are somewhat different in Tables [T] and from the rest 
of the paper. This is reflected in the sum of the absolute values of the weights, but has no effect on the 
various results other than scaling. 

First, we report on the algorithm LSQ. Each of the experiments in this case was repeated 30 times 
with data sets chosen randomly from the distribution fj,2 on S 2 . To test our algorithms, we computed the 
computed Gram matrix G COM given by 

The average maximum matrix norm of the difference between Q COM and the identity matrix of the 
same size indicates the error of the quadrature formulas. The results are shown in Table [TJ Based on 
these results we conjecture that in order to obtain stable quadrature formulas (i.e., with small condition 
number for the original Gram matrix Gn) exact for degree n > 1, one has to use at most 4d* +1 uniformly 
distributed points. In contrast, the theoretical guarantee in Theorem 14.1( c) requires 0(n 3q logn) points. 

As can be seen from the table, for a fixed degree n, the condition number k(Gjv) decreases as the 
number of points increases. For n > 140 and various sets of randomly generated points on the sphere, we 
do not obtain good numerical results. This might be due to a defect in the built in numerical procedures 
used by Matlab in computing the spherical harmonics of high degree at values close to —1 or 1. The 
situation was much better for the dyadic points; i.e., the centers of the dyadic triangles. 

For dyadic points on the sphere, the best result we obtained so far is n = 178 with 131, 072 points. 
As a further verification of this quadrature, we considered the following data. The data are constructed 
using coefficients {ae,k} for spherical polynomials up to degree 90, taken from model MF4 used for 
modelling the lithosphcric field. The model is computed by geophysicists at GeoForschungsZentrum 
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M 


n 


Error 


EKI 


min W£ 


max W£ 


pos 


k(G n ) 


A • 


Amax 


8192 


16 


2.41 * 10" 


15 


3.5449 


2.29* 10" 4 


7.88 * 10" 


4 


8192 


2.43 


0.607 


1.4730 




44 


4.32* 10- 


15 


3.5714 


-5.06* 10" 4 


0.0029 




8039 


37.52 


0.078 


2.8047 




64 


6.15 * 10" 


15 


5.5575 


-0.00664 


0.0073 




6068 


1695.1 


0.003 


3.9315 




84 


9.73 * 10" 


12 


82.152 


-0.16274 


0.1551 




4431 


3.52* 10 6 


2.6 *10" 6 


5.4851 


16384 


44 


4.43 * 10" 


15 


3.5449 


-9.50* lO" 13 


8.82 * 10" 


4 


16382 


9.19 


0.24036 


2.1590 




64 


5.25* 10" 


15 


3.5787 


-3.75* 10" 4 


0.0015 




16014 


51.9 


0.05964 


2.9150 




84 


7.10*10" 


15 


4.4757 


-0.0024 


0.0032 




13361 


944.86 


0.00612 


3.8457 




100 


1.94* 10" 


15 


9.1325 


-0.0077 


0.0072 




10625 


11896.1 


4.8* 10" 4 


4.6008 


32768 


44 


6.02* 10" 


15 


3.5449 


3.11 * lO- b 


2.90* 10" 


4 


32768 


4.270 


0.4157 


1.7652 




64 


7.09* 10" 


15 


3.5450 


-1.79* 10" 5 


5.23* 10" 


-4 


32761 


7.977 


0.8208 


5.2519 




84 


7.71*10" 


15 


3.5574 


-1.43* 10" 4 


7.92* 10" 


-4 


32410 


42.97 


0.0716 


2.7967 




100 


7.62* 10" 


15 


3.6777 


-4.28* 10" 4 


9.96* 10" 


-4 


30819 


145.6 


0.0250 


3.2967 



Table 1: The statistics for the experiments with the algorithm LSQ. M = \C\, n — 2 is the degree of 
spherical polynomials for which exact quadrature formulas were computed, N = n 2 , pos stands for the 
number of positive weights, k(Gn), A m i n , A max are the condition number, the maximum eigenvalue and 
the minimum eigenvalue of the matrix Gn respectively. 

Potsdam (Germany) based on CHAMP satellite data. We use those coefficients to construct the samples 
of a function / = ■ dt,kXt,k at the centers of 8 * 4 7 dyadic triangles. We then use our pre-computed 
quadrature based at these centers which can integrate spherical polynomials up to degree 178 to compute 
the Fourier coefficients ci^fe. The maximum difference between the vector {a^.fc} and the vector {a^fc} 
was found to be 6.66 * 10 -15 . 

Next, we considered the algorithm REC. In the context of spherical polynomials, the recurrence 
relations have to be chosen very carefully using the special function properties of the spherical harmonics 
Y^fc, in order to get stable results [55]. In the present situation, the polynomials tk have no special 
structure. Therefore, it turns out that the algorithm REC is not very stable for high degrees. However, 
when we took the centers of 8192 dyadic triangles as the quadrature nodes, and used the measure v TR as 
the starting measure, then we are able to obtain satisfactory quadrature formulas for degree 32. We note 
an interesting feature here that all the weights obtained by this algorithm are positive. These results are 
summarized in Table [5] below. 



n 


Error 


min(w^) 


max(w^) 




16 


4.196643* 10" 


14 


5.181468* 10 


-4 


2.538441 * 10" 


-3 


12.56637 


22 


5.302425* 10" 


13 


5.175583* 10 


-4 


2.543318* 10" 


-3 


12.56637 


32 


9.240386* 10" 


11 


5.154855* 10 


-4 


2.544376* 10" 


3 


12.56637 


42 


4.434868* 10" 


-8 


5.086157* 10 


-4 


2.544141* 10" 


-3 


12.56637 


44 


2.320896* 10" 


-5 


5.094771 * 10" 


-4 


2.562948* 10" 


-3 


12.56638 



Table 2: Quadrature constructed using REC on 8192 dyadic points 



Our second set of experiments demonstrates the local approximation properties of the operators 
<r n (C, W; h) for a smooth function h. For this purpose, we consider the following benchmark functions, 
considered by various authors [3ll [30l [T5J [10] , listed in (|5.1[) below. Using the notation x = (xi,X2,x 3 ), 
the functions are defined by 

, 9l (x) = ( Xl - 0.9)t /4 + (a* - 0.9) 3 / 4 , 

g 2 (x) = [0.01-{x 2 1 +x 2 2 + {x 3 -l) 2 )} + +cxp(x 1 +x 2 +x 3 ), 

,g 3 (x) = 1/(101-100^3), 

<? 4 (x) = l/(|asi| + l^al + lacal), 
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55 (x) 



cos 2 (^dist(x,(-l/2,-l/2,l/v5))), if dist(x,(-l/2,-l/2,l/V5)) < 1/3, 
0, if 



(5.1) 



In order to define the function h, we recall first that the B spline B 
m P- 131] by 



dist(x,(-l/2,-l/2,l/\/2)) > 1/3. 

of order m is defined recursively 



B m (x) 




Tfl — X 

-B m -x{x) H T B m -x(x - 1 

1 m — 1 



if m = 1, < x < 1, 
if m = 1, iel\ (0,1], 

if m > 1, a; e K. 



(5.2) 



The function B m is an m — 1 times interated integral of a function of bounded variation. We will choose 
h to be 



h m {x) = 2J B m (2mx-k), 



(5.3) 



for different values of m, in order to illustrate the effect of the smoothness of h m on the quality of local 
approximation. If m > 3, the function h m satisfies the conditions in Proposition 12 . 1 1 with S = m — 1. We 
note that the discretized Fourier projection operator aes(C, W; hi) has been called the hyperinterpolation 
operator [2"9"] . 

One example of the localization properties of our operators is given in the following table, where we 
show the error in approximation of gi on the whole sphere and on the cap K = §q 4510 ((— 1/V2, 0, — l/v^)- 
The operators were constructed using Driscoll-Healy quadrature formula [3] based on 4(n + l) 2 points, 
exact for integrating polynomials of degree 2n. The maximum error on the whole sphere, given in 
Columns 2 and 3, is estimated by the error at 10000 randomly chosen points; that on the cap, given 
in Columns 4 and 5, is estimated by the error at 1000 randomly chosen points on the cap. It is clear 
that even though the maximum error on the whole sphere is slightly better for the (discretized) Fourier 
projection than for our summability operator, the singularities of gi continue to dominate the error in the 
Fourier projection on a cap away from these singularities; the performance of our summability operator 
is far superior. 



n 


S2errhl 


S2errhb 


Kerrhl 


Kerrhh 


63 


0.0097 


0.0112 


3.4351 * 10~ 4 


6.5926* 10~ Y 


127 


0.0044 


0.0055 


8.0596* 10~ b 


6.5240* 10~* 


255 


0.0033 


0.0038 


1.4170* 10~ b 


1.1816* 10" 8 



Table 3: S2errhl — max x6 g2 |<?i(x) — cr n (C, W; hi, gi, x)|, S2errh5 = max xg § 2 |<?i(x) — 
cr„(C, W; /i 5 ,5i,x)|, Kerrhl = max xe K |fln(x) - <r„(C, W; hi, gi, x)|, Kerrhh = max x£S - |ffi(x) - 
<7„(C,W; /15, <7i,x)|, (C, W) are given by the Driscoll-Healy formulas. 



Theorem 13.21 points out another way to demonstrate the superior localization of our summability 
operator without an a priori knowledge of the locations of the singularities. Since each of the test 
functions is infinitely differentiable on large caps of different sizes, Theorem 13.21 suggests that the more 
localized the method, the greater is the probability that the approximation error would be smaller than a 
given number. To demonstrate also how our ability to construct quadrature formulas based on scattered 
data helps us to analyse the approximation properties of our summability operators, we took for the 
set C a randomly generated sample of 65536 points. For these points, the weights W computed by the 
algorithm LSQ yield a quadrature formula exact for integrating spherical polynomials of degree 126. 
We compare three approximation methods, the least square approximation from n| 3 , the approximation 
given by the operator <j§s(C, W; hi), and the approximation given by <J&s{C, W; hs). For each function, 
we computed the absolute value of the difference between the approximate value computed by each of 
the three methods and the true value of the function at 20, 000 randomly chosen points on the sphere. 
The percentage of points where the value of this difference is less than 10~ x is reported in Table [4] below, 
for x = 2 : 10. It is very obvious that cr 63 (C, W; h 5 ) gives a far superior performance than the other 
methods, due to its localization properties. 
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x — ► 


10 


9 


8 


7 


6 


5 


4 


3 


2 




51 








0.005 


0.02 


0.42 


4.44 


39.43 


94.79 


100 




LS 











0.04 


0.56 


5.32 


46.38 


95.45 


100 




55 


0.02 


0.19 


1.87 


16.89 


59.36 


68.34 


79.01 


93.09 


99.97 




51 





0.01 


0.09 


0.74 


7.94 


84.99 


99.19 


99.99 


100 


Go 


LS 





0.01 


0.15 


1.29 


13.29 


86.09 


99.28 


100 


100 




55 


0.39 


3.34 


41.95 


90.78 


94.52 


97.19 


99.18 


99.97 


100 




51 





0.01 


0.11 


1.26 


12.02 


91.87 


99.86 


100 


100 


Qi 
i/o 


LS 





0.01 


0.10 


1.49 


16.26 


93.12 


99.87 


100 


100 




55 


0.51 


5.43 


51.08 


82.22 


91.90 


95.79 


98.49 


99.87 


100 




51 











0.01 


0.18 


1.91 


18.28 


83.81 


99.97 


Qa 


LS 








0.01 


0.02 


0.25 


2.16 


21.24 


86.43 


99.98 




55 


0.01 


0.01 


0.04 


0.36 


3.47 


17.48 


40.98 


80.06 


99.88 




51 








0.01 


0.09 


1.12 


11.84 


88.94 


99.45 


100 


95 


LS 





0.01 


0.01 


0.15 


1.42 


15.23 


90.47 


99.75 


100 




55 


0.08 


0.64 


5.73 


66.82 


83.54 


88.74 


92.95 


96.78 


97.64 



Table 4: Percentages of error less than 10 _x for different functions, LS= Least square, Sl= error with 
063 (C, W; hi), S5= error with cr^C, W; h 5 ). For example, for the function g 3 , S5 was less than 10 - ' for 
82.22% of the 20000 randomly selected points, while SI (respectively, LS) was less than 10" 7 for 1.26% 
(respectively, 1.49%) points. 



x — ► 


5 


4 


3 


2 


3.0 


2.75 


2.5 


2.25 


51 


0.05 


0.635 


9.93 


97.45 





7.97 


92.75 


100.00 


LS 











100 





0.04 


30.93 


99.07 


55 


0.085 


1.015 


10.03 


97.49 


0.24 


51.97 


99.87 


100.00 



Table 5: Percentages of error less than 10~ x for e = 0.01, LS= Least square, Sl= error with cr^{C, W; hi), 
S5= error with aQs(C,W; h^). The random noise in the left half comes from the uniform distribution in 
[—e, e], that in the right half from the normal distribution with mean 0, standard deviation e. 

Next, we illustrate the stability of our operators under noise. Since our operators are linear operators, we 
assume for this part of the study that the target function / is the zero function contaminated either by 
uniform random noise in the range [— e, e], or a normally distributed random variable with mean and 
standard deviation e. We let C be a set of 65536 random points and computed corresponding weights W 
that integrate exactly polynomial up to degree 126. These were used in calculating a^{C, W; hi) and 
<763(C, W; /15) at each point of a test data set consisting of 20000 random samples from the distribution 
/i2- For each value of e = 0.1; 0.01; 0.001; 0.0001, the experiment is repeated 50 times and the errors are 
the averaged over the number of repetitions. The percentage of points at which the absolute computed 
value is less than \Q~ X is reported in Table[5]in the case when e = 0.01. The results for the other values of 
e were consistent with the linearity of the operator. We observe that in each case, both <jqz(C, W; hi) and 
<763(C, W; h 5 ) yield better results than the least squared approximation, while <7g 3 (C, W; h 5 ) is slightly 
superior to <Jsz{C, W; hi). 

Finally, we used our operator 022 (C, W; hi) with the MAGSAT data. Our purpose here is only to 
test how our methods work on a "real life" data. This data, supplied to us kindly by Dr. Thorsten 
Maier, measures the magnetic field of the earth in nT as a vector field. It was derived from vectorial 
MAGSAT morning data that has been processed by Nils Olsen of the Danish Space Research Institute. 
The measurements are averaged on a longitude-latitude grid with = 4° and A9 = 2° in geomagnetic 
coordinates. The radial variations of the MAGSAT satellite have been neglected in the dataset and, 
therefore, prior to the averaging process, the GSFC(12/83) reference potential model has been subtracted. 
The data results from one month of measurements, centered at March 21, 1980. We extract the East 
West component of the vectorial data as a scalar valued function on the sphere. Totally, there are 8190 
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data sites. A quadrature of degree 44 was computed based on those sites. Figure [2] shows the original 
data, its reconstruction using 022(C, W; hi), and the error in the approximation, \(T22(C, W; /17) — y\, as 
a map in the longitude-latitude plane. As can be seen from the figures, the reconstruction preserves the 
key features of the original data. 



6 Proofs 

In the interest of organization, we will prove the various new results in the paper in the following order. 
We will prove Theorem 14.21 first, since its proof does not require any preparation. We will then use 
Proposition 12.11 to prove Theorems 13.1 f a) and I3.2f a). Next, we will prove Lemma [6.11 and use it to 
prove parts (a) and (b) of Theorem 14.11 The remaining results in this paper involve probabilities. We 
prove Lemma 16 . 2 1 next . estimating the probability that the supremum norm of a sum of random spherical 
polynomials exceeds a given number. This lemma will be used immediately to prove Theorem I4.1f c). 
Finally, we will prove Theorems 13.1( b) , I3.2f b) . 

Proof of Theorem 14.21 It is convenient to prove (|4.9[) first. Since Pyu is a polynomial of degree 
Dk + 1, there exist real numbers ric,j(P) such that 

PVk= r k,j( P )Vj- 

Dj<D k + l 

Since the system {yu\ is orthonormal with respect to /1, 

r k,j(P) = / PykVjd/j,. 
Jn 
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If Dj < D k — 1, then the degree of Pyj is less than D k . Because of the lexicographic ordering where 
lower degree polynomials appear before the higher degree ones, this implies that Pyj £ V k -i- Since yk is 
orthogonal to Vk-%, it follows that rj t k(P) = if Dj < D k — 1. This completes the proof of (|4.9p . 

We observe that € span {iti, • • ■ , u p (fe)}. So, there exists a constant a such that ay p nA — Up(fc) S 
^(fc).!. Thus, af k y p (k) - /fcU p (fc) = af k y p (k) ~ u k+1 is linear combination of terms of the form f k Uj, 
1 < j < p(k) — 1- Since is the minimal index for which there exists a monomial with f k u p ( k ) € ^fc+ij 
each of the terms /fettj, 1 < j < p(fc) — 1 is in V4- It follows that otfkVptk) ~ u k+i € Vfe. Again, there 
exists a constant a' such that a'ufc+i — y k +i € Vfe. Therefore, writing f k — aa'fk, we conclude that 
hUp{k) - Vk+x = a'(afky p (k) ~ u k+i) + a'u k+1 - y k+1 e V k ; i.e., f k y p(k) = y k+ i - fkjVj- The 

Dj<D k 

first equation in (|4.8[) is now proved in view of (|4.9p . applied with p(k) in place of k, and the fact that 
-Dp(fc) = -Dfe+i — 1- We note that is a constant multiple of the monomial f k . The second equation in 
is proved in the same way. 
Using the second equation in (I4.8[) and (14. 9p . we obtain from the definition of c k ,j's that 

AkCu+ii — A k t k+ iy e d[i 



fktp(k)yedfi + y~] s ktj I t 3 yidn 



D k+1 -2<D 3 <D k 



V. 



= 51 n, m (f k ) / t p ( k )y m d^+ ^2 s k,jCj,e 

D t -l<D m <D t + l Jn D k+x -2<Dj<D k 

= ^ r e,m(fk)c p (k),m + 2^ s k,j c j,t- 

D e -l<D m <D e + l D k + 1 -2<D J <D k 

This proves (|4~T0|) . □ 

Next, we use Proposition 12 . 1 1 to prove Theorem 13. If a) and Theorem 13. 2\ a). 
Proof of Theorem I3.1f a). To prove part (a), we assume without loss of generality that n > 8, and 
let I > 1 be the largest integer with 2 £+2 < n. In view of (|3.4p , 

HZ-^^.W;^;/)^ < cE n/2 (f) < cE 2W (f) < c\\f - a^ih-J)^ 

OO 

< \\<T* 2k+1 (h;f)-a* 2k (h;f)\\ oc <c2- re <cn~ r . 

k=t+l 

This proves part (a). □ 

Proof of Theorem 13. 2[ a). Let K" be a spherical cap, concentric with K, K', and having radius equal 
to the average of the radii of K, K' . Let ip be fixed, C°° function that is equal to 1 on K" and equal 
to outside of K. Without loss of generality, we may assume that n > 8, and let I > 1 be the largest 
integer such that 2 i+2 < n. The direct theorem of approximation theory (cf. 26J) implies that there 
exists P £ IT^ such that 

M>--P||oo <c2- es . 

Therefore, using the definition of ||/||w r (i<")j we conclude that 



E 2W (M < \\fiP - Pa* 2i (h; f)^ < \\{f-a^(h;fM\ 00 + \\{il>-P)a^{h;f)\ 

"2.1 



< c{||/- ( t 2 %(/i;/)|| 0O ,k + 2-" S ||/|| 00 } 



< cl \\v* 2k+ i(h;f)-<7* 2k (h;f)\\oc,K + 2- nS \\f\\ 00 \ < cT 
Kk=i+1 



rt 



In view of (pO]) , 

||/-ff„(C,W;/i;/^)|| 00 , J r/ = \\M - <r„(C, W; h; /VOlUx' < \\M ~ <r n {C, W; h; /VOIU 

< cS n/2 (/V)<^ +1 (/^)<c2- rf <cn- r . (6.1) 
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E ^/(0(i-V>(0)<M/*;x-0 

feC\K" 



Since 1 - ?/>(C) = for C e K" , we may use (|2~8|) to deduce that for x £ K' , 

\a n (C,W;h;(l-i>)f,x)\ 

c(K,K',K") un ^ c(K,K',K") 
^ rs=5 11(1- WlU^kd < ^g=5 • 

Together with (|6.ip and the fact that r < S — q, this implies (|3.9|) . □ 

Next, we prove Lemma 16.11 describing certain extremal properties for the weig hts w^ SQ . These will 
be used in the proof of parts (a) and (b) of Theorem 14.11 



Lemma 6.1 Let n > 1 be an integer, N — <2£ , C be a finite set of points on S 9 and v be a measure 
supported on C. Let := v({(,}), £ £ C. 

(a) If the Gram matrix is positive definite, then the weights w^ S ® are solutions of the extremal problem 
to minimize ^2^c w c/ v i subject to the conditions that X)fec ^^(O = ^M' 

(b) If holds, there exist real numbers W^, £ £ C, such that \W^\ < for £ £ C and X^ec ^^P(0 = 
J Pdfi q for all P Gil'. 

Proof. In this proof, we will write G in place of Gat. The Lagrange multiplier method to solve the 
minimization problem sets up parameters Xe and minimizes 

Setting the gradient (with respect to w^) equal to 0, we get = v^J^e ^iVtiQ- Writing, in this proof 
only, Q = Xeyi, we see that — ffQ(£). Substituting back in the linear constraints, this reduces to 
Sjec v zQ(0yt{€) = These conditions determine Q uniquely; indeed, Q = ^ . G^yj. This proves 
part (a). 

The part (b) is proved essentially in j21j . but since it is not stated in this manner, we sketch a proof 
again. During this proof, different constants will retain their values. Let M = \C\, R M be equipped 
with the norm |||r||| = X^ec w d r d' ^ n this P r0 °f only, let S be the operator defined on II 9 by S(P) = 
(P(0)t£C € K M , and V be the range of S. The estimate 

[ \P\d t , q <c 1 J2v S \Pm (6.2) 

implies that the operator S : II 9 — * V is invcrtible. We may now define a linear functional on V by 



x*{v) = J S-\v)d^ q , re 



It is clear from (|6.2[) that the norm of x* is bounded above by c% . The Hahn-Banach theorem yields a 
norm preserving extension of this functional to the whole space R M . Identifying this functional with the 
vector (W^)^c, the extension property implies that X^{eC^^ > (^) = / Pd^q for all P £ LI', while the 
norm preservation property implies that \W^\ < c\v^ for (eC. □ 

PROOF of Theorem 14. II (a), (b). Let N = d^ +1 , r e R N , and P = YLi r tVl- In this proof only, we write 
G in place of Gn ■ Then 

r T Gr = n\ £ vmiOVmiO \ r m = ]T ^P(0 2 , 

and r T r = ||P|||. Therefore, (OJ) with p = 2 implies that c\r T r < r T Gr < c\v T v for all r £ R N . The 
statements about the eigenvalues of G are an immediate consequence of the Raleigh-Ritz theorem [T31 
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Theorem 4.2.2]. Using Lemma RTTT b). we obtain weights W$ such that X^ec W$yi(£) = ^M) ^ = 1> ' ' ' > 
and |W^| < cv^, £ £ C. During the remainder of this proof, we write u>£ = w^ S ® . In view of Lemma l6.1f a), 
we have 

£K| ^ \ £ w € 

This completes the proof of part (a). 

In order to prove part (b), we adopt the following notation during this proof only. Let I denote the 
N x N identity matrix. For any N x N matrix H, let \\H\\ denote sup ||i?r||, ||r|| — 1, r £ M N . We note 
that \\H\\ is the largest singular value of H . If H is a symmetric, positive definite matrix, then it is also 
the largest eigenvalue of H, and moreover, IrfiJr^l < ||i?||||ri||||r2||, ri,T2 £ ffi . Using (|4.4jl . it is easy 
to conclude using the Raleigh-Ritz theorem that ||G — I\\ < cn~ q , \\G~ 1 \\ < c, and hence, 




\G~ 



IG-^Z-G)!) < c|| C" 1 1| || G* — ^|| <cn- q . 



Let y(x) denote the vector (yi(x), • • • ,?/jv(x)) T for x S E> q . In view of the addition formula, ||y(x)|| 2 is 
independent of x, and hence, 

N 

• 1 



Consequently, we have 
u„ z " S( 9i 



< 
< 



y(x) T G- 1 y(e)^ 9 (x) 
J y(xf(G- 1 -/)y(e)^(x) 
y(x)||||y(e)||^ 



IG 



y(x) T y(£)^ g (x) 
1,0, •••,0) T y(£)| <cn- q n q +c<c. 



This completes the proof of part (b). 



□ 



The proof of the remaining new results in the paper are based on the following lemma, that gives a 
recipe for estimating the probabilities involving polynomial valued random variables. 

Lemma 6.2 Let n, M > 1 be integers, {ajj}jLi be independent random variables, and for j = 1, • • • , M, 
Zj = Z(u)j, o) g II* have mean equal to according to w,. Let B, R > 0, max |Zo (x)| < Rn q , and 

the sum of the variances of Zj be bounded by Bn q uniformly on S q . Lf A > and l2R 2 (A + q)n q logn < B 
then 

M 



Prob 



3=1 



> ^JY1B{A + q)n q logn < an' 



(6.3) 



Here, the positive constant c\ is independent of M and the distributions of ujj 



Proof. The proof depends upon Bennett's inequality [271 P- 192]. In this proof only, we adopt a slightly 
different meaning for the symbols L, V, r\. Let L,V,rj be positive numbers, and Xj, j = 1, • • • , M, be 
independent random variables. According to Bennett's inequality, if the mean of each Xj is 0, the range 
of each Xj is a subset of [-L, L], and V exceeds the sum of the variances of Xj, then for rj > 0, 



M 



Prob I £X,| > r] \ < 2exp 
3=1 



V 
1? 



9(Lr)/V) , 



(6.4) 
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where, in this proof only, g(t) := (1 + t) log(l + t) — t. We observe that g(t) = f Q J^(l + w) 1 dwdu. 
Therefore, if < t < 1/2, then for < w < u < t, (1 + w)" 1 > 2/3, and hence, g(t) > t 2 /3. Consequently, 
if Lr/ < V/2, then 

( M \ 

Prob I | Jfjl > r)\ < 2cxp(-77 2 /(3y)). (6.5) 

Now, let x £ S q . We apply (|6.5[) with Zj(x) in place of Xj, Rn q in place of L, Bn q in place of V, 
Tj = ^/3B(A + q)n q logn. Our condition on n ensures that Lrj/V < 1/2 with these choices. Therefore, 

M 



Prob | Z j( x )\ > ^/3B{A + q)n q logn < 2rT A ~ q . 



(6.6) 



Next, in the proof only, let P* = Y,f =l Zj, x* £ § q be chosen so that |P*(x*)| = ||P*||oo, C C § 9 be 
chosen so that |C| - cn q and <5 C (§ 9 ) < l/(2n). Then we may find f eC such that dist(x*,£* ) < l/(2rc). 
Since P* € n 9 t , its restriction to the great circle through x* and £* is a trigonometric polynomial of order 
at most n. In view of the Bernstein inequality for these polynomials [H Chapter 4, (1.1)], 



|P*(f) - P*(x*)| < nljPloodist^x*) < (l/2)|P*(x*)|. 



We deduce that 



Therefore, the event 



M 




M 


E^ 


< 2 max 
xec 


E^w 


i=i 


oo 





n is a subset of the union of the |C| events 



| ^(x)| > ^3B(A + q)ni logn, x e C. Hence, the estimate (jUJ) implies (jOJ) with d = 2c. □ 

We are now in a position to prove Theorem 14. 1( c). 

PROOF of Theorem 14. lf c). Let x £ § 9 . In this proof only, let = $4 n (fyx-£) — f 3>4„(/i;x-C)eZ/Xg(C). 
Then the mean of each is 0, and its variance can be estimated by 

Z?dn q (€) < f($4n(h; x ■ Ofd» q (0 < cn q . 



Finally, \Z^\ < cn q for each £. Hence, we may use Lemma [6.21 with cM in place of B, and c in place of 
P, to conclude that 



Prob 



sup 77 E* 4 "^*'^ _ / 



$4n(/i;X- CWq(Q 



> C 2 



?i9 log n 



71/ 



< cn 



and with M > cn q \ogn/r] 2 , 



Prob 



SU P T?E $ 4n(^;X-0- / *4n(ft;X-C)dA» a (C) > v] 



< cn 



Since any P £ H^ can be written in the form 



P(C) = / P(x)d> 4 „(fo;x-C)^ g (x), 
we see that with probability exceeding 1 — cnT A 1 



tec J 
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l ^(x)*4n(ft;C'X)d/* 9 (x) 

- J J P(x)$4n(/i; x • C)dM«(x)d/i 9 (C) 
< 77 / |P(x)|dp,(x). 



d/i g (x) 



For P € LI 9 , we may now apply this estimate with P 2 s n| n . □ 

Another immediate consequence of Lemma l6.2l is the following lemma, describing the the probabilistic 
behavior of the operator cr n (C, W; h). 

Lemma 6.3 Suppose that m > 1 is an integer, C — admits an M-Z quadrature of order m, and 

let W be the corresponding quadrature weights. Let R, V > 0, and for j = 1, ■ • • , M, ej be independent 
random variables with mean 0, variance not exceeding V , and range [— R, R]. Let g £ C{W), ||.<?||oo < 1, 
and E = {^j9(^j)}(jec ■ Then for integer n > 1 with (R 2 /V)(A + q)n q \ogn < c^m q , 



Prob \\a n (C, W;^^)^ > ci 



V(A + q)n q log) 
m q 



< C2U 



(6.7) 



Here, the positive constants C\, c 2 , C3 depend only on q but not on M and the distributions of e 7 



Proof. In this proof only, if £ = £j € C, we will write for ej and for the weight in W corresponding 
to £. We use Lemma [6.21 with = m q w^e^g(£,)^ n (h; £ ■ o), £ g C. We note that the random variable 
LOj in Lemma [6.21 is in this case. It is clear that the mean of each is 0. Since (|2.12[) implies that 
\w(] < cm~ q , (|2.7p shows that ||2?£||oo < cRn q . Moreover, for any x 6 S 9 , the variance of £^(x) does not 
exceed Vm 2q w^ n {h; £ • x) 2 . In view of the fact that are M-Z quadrature weights, (|2 . 1 2[) and (|2.T[) 
imply that 

^m 2 %;!$„0; £ • x) 2 < cm 9 ^ |u/ C |** (Ji; £ • x) 



< cm 5 



$ 2 (/i;C-x)d^(C) < cm q n q . 



Thus, we may choose B in Lemma l6.2l to be cVm q . The estimate ()6.7|) now follows as a simple consequence 
of ||B3J). □ 



We are now in a position to prove the probabilistic assertions of Theorems 13.11 and 13.21 
Proof of Theorem 13.1T b). To prove part (b), we use Lemma l(T3l with g = 1. Since the range of e/s 
is contained in [—1,1], we may take R = V = 1, and obtain from ()6.7() that 



Prob ||o- n (C,W;/i;E)|| 00 >C4 



(A + q)n q log n 
m q 



< c 2 n 



□ 



The choice of n with an appropriate C3, ensures that (A + q)n q logn/m 9 < n 2r . Therefore, 

Prob (\\a n (C,W;h;-E)\\oo > C ^ T ) < c * n ~ A - 

The estimate (|3 . 8() is now clear in view of (|3.4[) and the linearity of the operators cr n (C, W; h). 

Proof OfTheorem 13.2( b). We apply Lemma |6~B1 again with g = 1. As before, we may choose 
R = V = 1. The choice of n with an appropriate C3, ensures that (A + g)n 9 logn/m 9 < n~ 2r . Therefore, 
(|6.7p with these choices implies that 

Prob (IM^W^E)!!^ > c 4 n- r ) < c 2 n- A . 
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Together with (|3.9p and the linearity of the operators a n (C, W; h), this leads to (13.10|) . 



□ 



7 Conclusion 

We have described a construction of linear operators yielding spherical polynomial approximations based 
on scattered data on a Euclidean sphere. While the operators can be defined for arbitrary continuous 
functions on the sphere, without any a priori knowledge about the location and nature of its singularities, 
they are auto-adaptive in the sense that the approximation properties of these globally defined polyno- 
mials adapt themselves on the different parts of the sphere according to the smoothness of the target 
function on these parts. While the theoretical properties of these operators and their localization were 
studied in [TJ] , a bottleneck in their numerical construction was the construction of quadrature formulas 
based on scattered data, exact for integrating moderately high degree spherical polynomials. So far, it 
was possible only to compute quadrature formulas exact at most for degree 18 polynomials. We show 
that a simple-minded construction involving a Gram matrix is surprisingly well conditioned, and yields 
the necessary quadrature rules, up to degree 178. Using these newly constructed quadrature formulas, 
we are able to demonstrate that our constructions yield superior approximation properties to those of 
more traditional techniques of least squares and Fourier projection, in the sense that the presence of 
singularities in some parts of the sphere affects the degree of approximation by our operators on other 
parts far less than in the case of these other traditional techniques. We give probabilistic estimates on 
the local and global degrees of approximation by our operators in the presence of noise, and demonstrate 
its use in the modeling of a "real life" data set. We also describe a theoretical algorithm for construction 
of data dependent multivariate orthogonal polynomials and their use in the construction of quadrature 
formulas, analogous to the univariate algorithms in the book [TT] of Gautschi. 
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